Done.
[ NOTE, FIGURE ]: Red horizontal line indicates a signed R^2 of 0.9
Setting soft-thresholding power to: 8.
Power-transforming the gene-gene similarity matrix...Done.
---------------------------------------------------
4. Convert into topological overlap matrix (dissTOM)
---------------------------------------------------
Creating dissTOM...Done.
Performing hierarchical clustering on dissTOM...Done.
---------------------------------------------------
5. Identify modules (clusters)
---------------------------------------------------
Merging modules that have a correlation ≥ 0.9 ...Done.
[ NOTE, FIGURE ] Plotting identified clusters before and after merging.
Visualizing a simplified representation of the circadian GCN
Annotate the network
Identify rhythmic modules
Show code
# Obtain a list of genes in each modulel_module_genes <- mods[["module_genes"]] |>arrange(module_identity) |>group_split(module_identity) |> purrr::map(~ .x |>pull(gene_name) ) |>setNames(unique(mods[["module_genes"]]$module_identity))# Load results of rhythmicity analysesdb_rhy <- purrr::map( sample.names,~load_rhy_genes(sample = .x )) |>setNames(sample.names)# Set your p-value of choice###-###-###-###-###-###-###-###-# col_pval = "BH.Q"col_pval ="default.pvalue"# col_pval = "raw.pvalue"###-###-###-###-###-###-###-###-# Obtain a list of rhythmic genes in each tissuel_rhy_genes <- purrr::map( sample.names,~ db_rhy[[.x]] |> purrr::map(~ .x |>filter(if_all(all_of(col_pval),~ .x <0.05 ) ) |>filter( ID %in%unlist(l_module_genes) ) |>pull(1) |>unique() ) |> purrr::compact()) |>setNames(sample.names)
Modules vs. Rhythmic genes
Show code
# LIST ONE - WGCNA moduleslist1 <- l_module_genessapply(list1, length) |>print()
“We begin by calculating the intramodular connectivity for each gene. (In network literature, connectivity is often referred to as”degree”.) The function intramodularConnectivity computes:
the whole network connectivity kTotal,
the within (intra)module connectivity kWithin,
the extra-modular connectivity kOut=kTotal-kWithin, and
the difference of the intra- and extra-modular connectivities kDiff = kIn - kOut = 2*kIN-kTotal
Show code
# From what I can tell, colorh1 in the tutorial refers to moduleColorscolorh1 <- mods[["modules"]]$colorsadj_matrix <- mods[["adj_matrix"]]# Calculate the connectivities of the genesAlldegrees1 = WGCNA::intramodularConnectivity(adjMat = adj_matrix,colors = colorh1) |>mutate(gene_name =rownames(adj_matrix),across(matches("^k"),~round(.x, 2) ) ) |>glimpse()
# tidydata.db |> # purrr::map_int(# ~ nrow(.x)# )# We work with two sets:nSets =2;# Find the common set of genes across all samplescommon_genes <- tidydata.db |> purrr::map(~ .x |>pull(gene_name) |>unique() ) |> purrr::reduce( intersect ) |>intersect( mods[["module_genes"]]$gene_name ) # filter to keep only genes used in creating REF networksamples.list <- tidydata.db |> purrr::map(~ .x |>filter( gene_name %in% common_genes ) )test_samples <-setdiff(sample.names, which.sample)l_mp <- purrr::map( test_samples,function (x) { multiExpr <-prep_data_module_preservation(data = samples.list,ref = which.sample, test = x )# CHECKS#---#---#---#---#---#---#---#---#---#---#---#---# Check that the data has the correct format for many functions # operating on multiple sets: exprSize = WGCNA::checkSets(multiExpr)# Check that all genes and samples have sufficiently low numbers of # missing values. gsg = WGCNA::goodSamplesGenesMS(multiExpr, verbose =1);cat("All samples okay?\n", gsg$allOK)#---#---#---#---#---#---#---#---#---#---#---#---# prep data#---#---#---#---#---#---#---#---#---#---#---#---# Expression data multiExpr_1 =list(ref =list(data = multiExpr[[1]]$data), test =list(data = multiExpr[[2]]$data) )## filter to keep only the genes that we are working with mod.identity <- mods[["module_genes"]] |>select( gene_name, old_labels ) |>filter( gene_name %in% common_genes ) |># !!this step is necessary!! #arrange(gene_name)## specify the module identity of the genes moduleColors <- mod.identity |>pull(old_labels) multiColor =list(ref = moduleColors)# Run module preservation#---#---#---#---#---#---#---#---#---#---#---#--- mp = WGCNA::modulePreservation( multiExpr_1, multiColor,referenceNetworks =1,nPermutations =200,calculateQvalue =FALSE,quickCor =0,verbose =1 ) mp })
Flagging genes and samples with too many missing values...
..step 1
All samples okay?
TRUE ..checking data for excessive amounts of missing data..
..calculating observed preservation values
..calculating permutation Z scores
..Working with set 1 as reference set
Flagging genes and samples with too many missing values...
..step 1
All samples okay?
TRUE ..checking data for excessive amounts of missing data..
..calculating observed preservation values
..calculating permutation Z scores
..Working with set 1 as reference set
Flagging genes and samples with too many missing values...
..step 1
All samples okay?
TRUE ..checking data for excessive amounts of missing data..
..calculating observed preservation values
..calculating permutation Z scores
..Working with set 1 as reference set